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Abstract 



Important gaps remain in our understanding of the thermodynamics 
and statistical physics of self-gravitating systems. Using mean field the- 
' O ■ ory, here we investigate the equilibrium properties of several spherically 

G ' symmetric model systems confined in a finite domain consisting of either 

j^ , point masses, or rotating mass shells of different dimension. We establish 

a direct connection between the spherically symmetric equilibrium states 
of a self-gravitating point mass system and a shell model of dimension 
3. We construct the equilibrium density functions by maximizing the en- 
^ ' tropy subject to the usual constraints of normalization and energy, but 

in I we also take into account the constraint on the sum of the squares of the 

individual angular momenta, which is also an integral of motion for these 
symmetric systems. Two new statistical ensembles are introduced which 
incorporate the additional constraint. They are used to investigate the 



(<-^ ■ possible occurrence of a phase transition as the defining parameters for 

(^ ' each ensemble are altered. 

1 Introduction 

'■^ ' The observation that a number of different types of astronomical objects appear 

W . to be in thermodynamically relaxed states has motivated theorists to under- 

stand the thermodynamics and statistical physics of self-gravitating systems. 
Of particular note are the globular clusters, consisting of about a million stars. 
Besides having relaxed cores, these structures appear to be organized in two 
k^ , distinct classes characterized by radically different density profiles, x-ray pro- 

^ • duction, and other features. Q] This suggests that globular clusters may exist 

. 5t 1 in different thermodynamic phases. However, in contrast with normal "chemi- 

cal" systems, which have been successfully described by thermodynamics at the 
macroscopic level, both the infinite range and short distance singularity of the 
Newtonian gravitational potential introduce problems in the statistical theory 



of phase transitions which make their analysis a challenging task. The descrip- 
tion of the system can be simplified by going to the Vlasov limit, i.e. by letting 
the number of particles become large while controlling the total mass, M, and 
energy, E. In this limit the system is described by the single particle density 
f(x,v,t) in the /z (position, velocity) space, which is employed by most of the 
standard treatments, including the present work.[E| We refer to this reduced de- 
scription as mean field theory (MFT). While MFT avoids the problem of dealing 
with an N body formulation, the difhcultics introduced by the singularity and 
long range of the potential persist. 

In the early sixties Antonov investigated the equilibrium behavior of iso- 
lated gravitational systems in MFT.[|| To circumvent the problem of escape, 
he confined the mass to a finite region by introducing a rigid wall. By fixing 
the total mass and energy, he showed that maximum entropy solutions for / 
are spherically symmetric in position and have the expected Maxwellian ve- 
locity dependence. However, he proved that there is no global maximum to 
the entropy: while extremal solutions can exist, these are at best local max- 
ima. Antonov |^ , as well as Lynden-Bell and Wood [Q , closely investigated the 
spherically symmetric system confined in a sphere and showed that when the 
mass and energy are fixed there are no entropy extrema above a critical radius 
(i? > — 0.335GM^/£') . When the radius is less than this value, the stability 
of the extremum solutions were studied by several authors (Lynden-Bell and 
Katz §, Katz § §, Padmanabhan g §, and Bavaud |l§). They found that, 
in general, above a critical density contrast [p (0) / p [R) > 709) all extrema 
are unstable, i.e. they are not local maxima. Lynden-Bell and Wood coined 
this phenomena the gravothermal catastrophe and it is also referred to as the 
Antonov instability. In such a system, there is no upper bound on the entropy 
and a state of arbitrarily large entropy can be constructed from a centrally 
concentrated density profile by shifting more of the mass towards the center 
(core- halo structures have higher entropy). 

More recently, Kiessling has investigated the thermodynamic stability of the 
full N body point mass system confined in a spherical box using the canonical 
ensemble. To avoid the short range singularity he regularized the Newtonian 
interaction by softening it (letting the potential approach a finite value as the 
origin is closely approached). He showed that in the limit that the softening 
vanishes (Newtonian interaction), the canonical equilibrium measure is the su- 
perposition of Dirac measures at any temperature, meaning that the system is 
in a collapsed point mass state. We have to emphasize that this is the equi- 
librium solution when the system is in thermal equilibrium with a heat bath. 
Also, based on the results of the finite N-particle system, with proper scaling 
of the particle mass as we take the mean-field limit, he showed that the single- 
particle density function is proportional to the Dirac distribution. Therefore 
the system's equilibrium state is the collapsed state in the canonical ensem- 
ble. Kiessling's conclusions do not contradict the earlier work applying MFT to 
the microcanonical ensemble (fixed mass and energy) described above, since no 
global entropy maximum was found in that case either. 

In fact, the pure Newtonian potential is never a correct picture in general 



because of the finite size of stars and atoms. In a nearly hidden appendix of 
their paper describing the gravothermal catastrophe, it was first pointed out by 
Lynden-Bell and Wood that if we modify the singular \/r Newtonian gravita- 
tional potential at the center by introducing a small minimal distance between 
the particles (the so called hard-sphere model) that complete collapse would be 
avoided and a global entropy maximum should exist. [0 They further conjec- 
tured that in this situation a first order phase transition to a centrally concen- 
trated core-halo configuration would occur as the system energy is reduced. Sev- 
eral authors demonstrated the existence of a phase-transition in special mean- 
field models with a modified gravitational potential. Hertel and Thirring |Q 
were the first to show analytically that a gravitating system can undergo a first- 
order phase transition. Although their system is not purely gravitational and 
has no singularity, the pair-interaction potential is purely attractive and has a 
fixed value when the pair of particles are in a given sub-domain. Also Lynden- 
Bell and Lynden-Bell ||l^ showed the occurrence of a first-order phase transition 
in their special gravitational system, consisting of point particles distributed on 
a shell which cannot shrink into a point mass (inner boundary) or expand to 
infinity (outer boundary). Kiessling et. al. also applied the hard-sphere model 
[O to study planet formation: They were able to explain the existence of ob- 
servable planets by showing that the mass belonging to the condensed phase is 
well below the Jeans mass.pl For finite N-body systems, a gravitational first- 
order phase transition was first observed dynamically by Youngkins and Miller 
[ISl. They investigated a model consisting of irrotational, concentric, spherical 
mass shells confined between two rigid spherical boundaries. The system was 
studied both theoretically in the mean-field limit and numerically by N-body 
simulations in the microcanonical, canonical, and grand canonical ensembles. 
The analysis for this one dimensional system showed that the system under- 
goes a first-order phase transition instead of a gravothermal catastrophe. As 
expected in gravitational systems, there were some discrepancies in the results 
for different ensembles; however the numerical N-shell simulations were always 
in good agreement with the corresponding mean-field predictions. 

Much earlier, Eddington ||l^ determined the form of the general station- 
ary solution of the Vlasov equation for a spherical system with an anisotropic 
velocity distribution which obeys the Schwarzschild law [Q . This model explic- 
itly depends on the square of the angular momenta P , but it also includes the 
isothermal part: /(r,v) a e^'^^e"'*'' . The model was presented much earlier, 
but it may have been forgotten. Later, several phenomenological models were 
improved by including Eddington's anisotropic term (King-Michie models and 
others ||lj [00), giving a better fit to the observed density profiles of globular 
clusters. However, in some cases, a good fit was not obtained for globular clus- 
ters with core-halo structures. In other words, a fair number of globular clusters 
which have a small and very dense core surrounded by a thin halo structure, do 
not obey these empirical density-fit models. On the other hand, in addition to 
the system energy, the sum of the angular momenta squares L2 = X) 'f i^ ^^^^ 
an integral of the motion for any isolated spherically symmetric system in the 
mean- field limit B], and should not be ignored. 



The purpose of our work is to present a more general approach for investi- 
gating the equilibrium properties of confined, spherical, systems then those of 
Antonov, Lynden-Bell and Wood, etc. mentioned above which takes into ac- 
count both integrals of a spherical gravitating system, and to introduce idealized 
dynamical "shell" models which also satisfy these constraints. Model gravitat- 
ing systems consisting of a collection of concentric, infinitesimally thin, spherical 
shells were first introduced by Henon pq | . They are useful for investigating the 
initial stages of evolution of a spherically symmetric self gravitating system, be- 
fore the onset of binary formation arising from three body effects. |l| They have 
the further advantages of ease and accuracy of algorithm construction, since it 
is possible to analytically solve for the motion of each shell between encounters, 
eliminating the need for the tedious and slow, step-wise integration of coupled, 
nonhnear, differential equations. Q ||l^ pO|| . 

In the present work we consider the mean field theory of a system of gravi- 
tating point particles moving in three dimensional space, as well as that of thin, 
rotating, spherical mass shells with angular momentum vectors restricted to 
manifolds of one, two and three dimensions. We first determine conditions for 
the equilibrium one particle probability density function f {x,p)oi a unit mass 
particle (shell) by finding the entropy extrema with respect to the constraints 
of (1) the normahzation, (2) system energy E, and (3) the sum of the squares 
of the angular momentum, L2. We then show that the introduction of the new 
integral L2 suggests a new type of canonical ensemble (T — 7), in addition to 
the extension of the microcanonical ensemble {E — L2) ■ A nonlinear differential 
equation governing the radial density valid for each ensemble is derived for the 
case of the three dimensional point mass system, and each shell system. We 
then prove that the radial density of the shell system with angular momentum 
confined to the Euclidean plane satisfies the identical differential equation as the 
three dimensional point mass system, and we carefully study the equilibrium 
solutions for this case numerically. The stability of the extremum solutions 
of each model is investigated in both microcanonical {E — L2) and canonical 
(T — 7) ensembles. At first glance it is natural to anticipate that the centrifugal 
barrier associated with the additional constraint will eliminate any tendency for 
complete core collapse without introducing an inner boundary in the system 
or changing the gravitational potential by other means. We conclude by in- 
vestigating the possible presence of a phase-transition which would remove the 
gravothermal catastrophe . 

2 The Entropy Extrema 

2.1 Spherically symmetric point mass systems 

The primary goal is to evaluate the equilibrium single particle probability den- 
sity function / (a;,p) which maximizes the entropy in the mean-field limit. Con- 
sider the spherically symmetric isolated point mass system in three space di- 
mensions with total mass M confined in a sphere of radius h. We choose units 



where G — M = b = I and introduce spherical coordinates x — (r, ip, d) . The 
Lagrangian per unit mass of a single particle moving in the mean field potential 

Hr) (D is 



L 



Y + V ( 



d +Cp^ sin^ ?9 - $ (r) , 



(1) 



where, for a Newtonian pair-wise interaction, $ (r) is given by 



^ir) = - f f G(r,r/) +G(r/,r) f U , p'\ cP x^ dP p' , 



(2) 



where 



G 



[r, r^j 



A ejr-rQ 



and, as usual, 6 (r) denotes the Heavyside step function. The Hamiltonian is 
then 



H = lv^ 



l^, 



/2 



2r2 2r2 sin^ d 



+ $(r) 



(3) 



where p = {v,lip,ls) are the corresponding canonical momenta. Our plan is first 
to construct the entropy extrema, and then verify whether or not the solutions 
are local maxima. The entropy of the system isyjg 



S = - 



JlnffFxcPp. 



(4) 



and the constraints for which we need to find the extremum are: (1) normaliza- 
tion of /, (2) energy conservation in the complete system, and (3) conservation 
of L2 : 



1 = 



fcPxcTp 



(5) 



where 



E 



ff(fx(fp 
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f[-v' + ^l' + -^ir)]d'xd'p 



(6) 
(7) 



l' = 



li + -^ 



sin 1} 



Introducing Lagrange multipliers a, f3, 7 we have an extrcmuni for the functional 
iS* when 

= 6{S-a~PE--/L2). (8) 

Taking the first variation of each term in (||), and asserting (||), (Q), (H), and 
(H), we obtain: 



Iv^ + ^f + Hr) 



+ -fl' } Sfd-'xd-'p 



^-.2 , ^f„\\ ( P , _, \ ;2 



0--y j |(ln/ + l) + a + /? 
from which finally 

= - (In/ + f ) - a - /3 l^v' + $ (r) 1 - I ^ + 7 1 P. (9) 

Thus the one-particle probability density function (pdf ) is 

/ = exp[- (a + 1) - /3 Qz;2 + ci>^ _ (^A + ^^ ^2]^^ (^q) 

In order to obtain the radial density p (r) , we have to integrate / over the 
other variables. To ensure that the integrals over v , l^, and l^p converge, the 
following conditions are necessary: /? > and /3/2r^ +7 > at any r. Therefore 
the second condition is 2^/P > —1/6^ = —1. Using K = exp [— {a + 1)] , we 
get 

p{r)= J J fd^pdipd^ = J J /^y^^^exp(-/3<l>)dv5d^ 

and, with the Poisson equation in a spherical coordinate system, 

where p is the volumetric mass density. In some cases, it's more convenient to 
use the linear (radial) density instead of the volume density: 

Note that because M — \, the radial probability density function and the linear 
mass density function are the same. Introducing a new function ^ = /3$ and 
employing |ll| we can rewrite (|l2| ) as 



i("'?)-«^-<-)'v^(^-^)"' 



e-*M 



obtaining a closed equation for ^. This in turn can be simplified by introducing 
constants C and F, 



C = K{2^)i^>Q, (14) 

27 

This is the final form of the differential equation for the scaled potential which 
we will solve using numerical methods. Finally, we can evaluate E, L2, and S 
in terms of the Lagrange multipliers, the density p (r), and the potential $ (r) . 
After integrating (|^), (^), and (|^), we obtain: 

(15) 
(16) 



s 


5 

= "+2 


+ /5 / p^dr 
Jo 


L2-- 


.6 / 
= / P 

Jo \ 




-1 
dr 


1 


Jo 


^ 1 


+ 2 


w 


V^ + 27r2 



E=W-a+ P\ IT-KZT^ + o M^ (17) 



2.2 Shell systems 

In this section we consider the system of spherically symmetric, infinitesimally 
thin, mass shells confined in a sphere with radius h. We can define three basic 
types of the model with dimensions d = 2, 3, 4. The one-dimensional d=\ non- 
rotating shell system is discussed in ||T^ . In the d = 2 case, every shell rotates 
about a fixed axis. When d = 3, the rotational axis of every shell is in a fixed 
plane, while d = 4 means that every shell can rotate about any arbitrary axis. 
Therefore we use d — 1 angles as coordinates. Let us consider these systems 
in the mean- field limit again using units where M — G = h = 1. With the 
potential discussed above (0), the Lagrangian and the Hamiltonian of a unit 
mass shell in a spherical coordinate system are, respectively,: 

1 1 " 

^=-/ + -/Y.^l-^^) (18) 

/c=l 



1 3V" P 

4" 



H=-v^ ^- ^^=}- ^ + $ (r) (19) 



where (pj. {k — 1,2,3) are the angles around x, y and z, the Ik are the x,y,z 
components of the angular momentum per unit mass, and we have used the fact 
that the moment of inertia of a shell with unit mass is 2/3r^. In the equations 
above, we use n = d — 1 which is the number of degrees of freedom coming 
only from rotation. This is a generalization of the previous model, and so is 
the method to find the equilibrium pdf, f {x,p). As in the previous section, in 
order to get the equilibrium solutions first we have to find the cxtrcmum of the 
entropy: 



S 



f In fd"+'xd"+'p 



with respect to the three constraints of normalization, E, and L2 '■ 

1 = / / fd^'+^xd'^+^p, 



(20) 



Lo = 



f(J2ll]d-+'xd-+' 



P, 



\k=l 



(21) 



E 



f I i„2 ^ 3 ELi^fc ^ 1 ^ (^) ) d^+^xd^+^p. 
2 4 r^ 2 



(22) 



The solution for the variational problem can be obtained easily, and the one- 
particle density function is now 



/3( -1-2 + $ 



cxp 



-I^+"']ll^ 



f = exp [— (a + 1)] exp 
Therefore the radial mass density function is 
p{r) 

From the Poisson equation, again using using ^ = /?$, we find 

3/3 



2 

k 
fe=l J 



(23) 



/d"+ V> = K (27r) '^^(^^+^] exp (-/?$) . (24) 



^rf)-^(2'^)'^^'^(S+^) "^^p(-*) (2^) 



which can be simplified to 



d_ / ,d*\ 
dr \ dr J 



2"-)=cf4T + r) e-*M 



(26) 



C = —L^TT^^^K^/P > 0, 



im- 



(27) 



% > ' (^«) 



Comparing the results (26) for the 3 dimensional shell system (n = 2) to that 
of a 3 dimensional point mass syste m (|l^ ), we can see that these two systems 
are equivalent in so far as ( [l3| ) and (^6[) have the same form. 

It is useful to evaluate S,L2, and E in terms of the Lagrange multipliers, 
p{r), and $ (r) . The integration of @, (|l|), and (||) yields 

S^a + -^ + 13 / p^dr (29) 



f/'d^-f* (30) 



1 f^ fin 1 $ 



E= — + p[—^ + -]dr. (31) 

3 T-7 ensemble 

In the previous sections we derived expressions for the entropy extremum solu- 
tions for our model systems in terms of the local radial density. A third Lagrange 
multiplier,7, was introduced to satisfy the new constraint on L2 . Thus speci- 
fying E and L2 (and, of course, M) defines the analogue to the microcanonical 
ensemble for these systems. Alternatively, by fixing /3 and 7, we can define the 
analogue of the canonical ensemble, in which E and L2 are not fixed, but their 
average is determined by /3 and 7 , where 7 = ^-. We call this the T — 7 
ensemble. It can be modeled by imagining that the system is in contact with a 
heat bath with constant /3 and also an P bath at constant 7. The new P bath 
corresponds to the fact that we allow some P exchange between the system 
and the bath. Since the system is spherically symmetric in position, while the 
system can also exchange angular momentum with the bath, its vector average 
will vanish. As an example, we can imagine a globular cluster which is embed- 
ded in some large spherically symmetric stellar neighborhood with an isotropic 
velocity distribution. The mean angular momentum of both system and bath 
is zero, and only P and energy exchange can occur (for the moment, we do not 
take into account the possibility that particles can escape from the cluster). 

In performing calculations it is more convenient to use the T — 7 ensemble 
than the E — L2 microcanonical ensemble because in the latter we have to find 
the Lagrange multipliers from the given E and L2- As can be seen from ^, 
|l6| , and |l7| this is a nontrivial and laborious task. In this ensemble the relevant 
thermodynamic potential is an extension of the Helmholz free energy: 

F = E-i^S+lL2 (32) 

and equilibrium states, if they exist, minimize F. 



4 Stability 

From the variational problem, we only know the extremum solutions. In order 
to separate the unstable solutions from the locally stable, we use the modified 
method of Poincare's linear series of equilibria. Following Katz, QQ] we can 
generalize the method to the case of a functional. From now on, we discuss 
the generalized version of the method which has to be applied to our models to 
determine whether an extremum solution is stable or unstable. We outline the 
method below. Let's assume that we want to find the maximum of the functional 
F* (/, s) which depends parametrically on s. The function / : f2 — > i?, where 
f2 is a compact domain 17 e i?'^, and the local maximum of the functional is 
the stable solution. Suppose we partition il and consider the vector x e i?" 
with elements xi = f [yi) (yi £ fl is not a coordinate but an indexed element 
in Q). Instead of dealing with the functional F*, we can construct a function 
F : i?" -^ R such that F {x, s) « F* (/, s). The accuracy can be controlled by 
refining the partition (increasing n). The problem is then shifted to finding the 
extrema of F: 

d,F {x, s) = (33) 

Denote the extremum solutions of the problem by a; = {Xa (s)} where a,b = 
1...N, labels different extremal solutions. Assume that the first and second 
derivatives of F are continuous in x and s, and Xa is continuous as well. As- 
sume also that the matrix {—didjF) has a non-degenerate eigenvalue spectrum, 
which we may consider to be ordered: fcia (s) < fc2a (s) < • • •< kna (s) < ..., 
and further that {—didjF) is diagonal. (We can always transform it into that 
form.). Let's evaluate F at the extremum points x = Xa (s) where 

d,F{Xa,s)^0. (34) 

Then, as the parameter s is varied, on the extremum labeled a, 

= -d,F {Xa, s) = (5,5,F)„ + J2 {dAF)a ^i = {dsd^F)^ - ha (s) K, 



(35) 






The stability of the extremum is determined by the second derivative of -Fa (s) — 

F{Xa,s), 

F^ - id'.F)^ + Y: idAF)a K - {d^F)^ + Y: ^^^. (37) 



1=1 



The stability will change only when one of the kia changes sign, and a change 
in stability occurs only at bifurcation or limit points p] M . Therefore we have 



10 



to investigate the dependence of Fa on s in order to decide how the stabiHty 
changes at an a — 6 bifurcation or hmit point . If we look at (B7h , one can see 
that when a particular kia (sq) changes sign from positive to negative, at that 
point lim Fa — +00 and similarly, from the other b branch, lim Ft — — cx). 

Going back to our original problem, we have to apply the method to our 
model systems both in the E — L2 and T — 7 ensembles. The only difficulty 
is that the method discussed above only allows us to include systems with one 
control parameter, s. In our case, we are free to fix either E or L2 in the E — L2 
ensemble and /? or 7 in the T — ^ ensemble and regard the other as the control 
parameter. Afterwards, we can change the previously fixed parameter and apply 
the method again for a wide range of parameter sets. In the E — L2 ensemble a 
natural choice of parameter is E with one fixed value of L2 , as the entropy has 
to be a local maxima if the system is in a locally stable state. If we are at the 
extremum solution points 

For a different value of L2, we can clearly see which branch of the extremum 
solutions are unstable. To obtain a complete description we can use the same 
method if E is fixed: 

A dS dS 
dL2 dL2 

In the T — 7 ensemble, first consider the case where we fix 7, and we are looking 
for the maximum of the functional —(3F . On a branch of the extremum solutions 

dF _dE ^ S I dSdE I dS dL2 7^,7 dL2 



d(3 d(3 /3^ (3 dE d(3 (3 dL2 d(3 /3^ ^ f3 d(3 
■m^P and ^ 



Using ^ — 13 and -^ = 7, we easily find 



dF _ S 7 

d{(3F) ( dF\ , , 

' F + /3— - = -E. (41) 



df3 \ ' d{3 

If we construct the stable branches for several values of 7, we can build up a 
general picture of the stability of the extremum solutions. We can apply the 
method for the case of a fixed /3 as well. The extremum solutions are stable 
when —13F is maximum: 



dF 

d7 


dE 
97 


IdSdE 1 dS dL2 I 7 9^2 
f3dE d-f f3dL2 d-f ^ P ' ^ [3 dj 

dF 1 
dl - /^' 
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From the results derived above we have an easy to apply tool to separate 
the unstable solutions. In the microcanonical {E — L2) description, we have 
to inspect the extremum solutions in the j3 — E plane for several fixed values 
of L2 or in the 7 — L2 plane for several fixed values of E. However, in the 
T — 7 ensemble, we have to inspect the extremum solutions in the {—E) — [3 
plane for fixed 7 values, or we have to consider the extremum solutions in the 
(— L2) — 7 plane for the case of fixed (3. Either way, as we will show below, we 
can generalize the approach to two parameters. 

5 Numerical method 

In order to find the entropy extremum solutions for the systems above, we have 
to solve ( [l3| ) or (|2g), for the 3D point mass system, or the three types of shell 
system. Generally, all of the differential equations can be written in the form of 



dy 
dr 




d^ 


y 


dr 


J.2 



-■^(r) 



(44) 



where ^ = /?$ and n = d — I (as above, d is the dimension of the system). 
Of course, the parameters F > — 1 and C > are different case by case. But 
it's quite interesting to mention that these systems behave similarly, regardless 
of whether we deal with a point mass system or a shell system, as long as the 
dimensions are the same. This is no longer a surprise if we look back at the 
Hamiltonians. But the similarity does also mean that, for example, a three 
dimensional point mass system is equivalent to a three dimensional shell system 
as far as the one particle density function and the stability of solutions are 
concerned. Also another advantage of the analogy is that we can dynamically 
model a three dimensional point mass system with a 3 dimensional shell system. 
Both systems should show the same equilibrium properties in the mean-field 
limit. Of course, we have to reassign the momentum of inertia of a shell to a 
different value in order to get exactly the same Hamiltonian in both cases. 
Equation (B) should be solved with the following boundary conditions: 

y{0) = (45) 

*(1) - -/? 

Unfortunately, there are two things which make finding the solutions more diffi- 
cult. First of all, (^ ) has a singularity at r = and, secondly, the existence and 
uniqueness of the solution are questionable for any given C, /3, and F. We cannot 
simply set C and F to satisfy our constraints of specified E and L2 because we 
don't have explicit forms of E and L2 in terms of C and F : the constraints are 
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functionals of p and $. To eliminate this problem we can change our boundary 
condition problem to an initial value problem because, in the latter case, we can 
ensure the existence and uniqueness of our solution for any C > and F > — 1. 
Therefore we choose 

y(0) = (46) 

*(0) - *o 

where ^o G {^oo, +oo), and we are able to construct the solution around r = 
in the form of a power series of $ and y. From the numerical point of view, we 
should take the power series of ^ and y in ri (a sufficiently small radius around 
r = 0), and then continue the integration of ( [44|) numerically from that point 
with the new initial values 

y{ri) -= 2/01 (47) 

*(ri) ^ *oi 

Although the solution we get from (Efl) is unique, it does not satisfy both (Bsh 
and normalization. However, we can prove that we can't find any Xi,X2 G ^o 
such that 'I'l = ^2 + c, where Vq consists of those solutions x — (2/i *) which 
satisfy the same initial value problem of (E3) and (Mfl) with C = Co > and 
any F > — 1. Therefore all of the solutions are unique with fixed C and F or, in 
other words, all x G Vb are unique. Now let's consider Vi, which can be defined 
the same way as we defined Vb, but with C = Ci ^ Cq. We can also prove that 
for any X2 ^ ^i' there is only one Xi ^^0 '■ X2 — Xi~ (0; c). Thus we can find 
all of the physically unique solutions by solving the initial value problem with 
a fixed Cq. In practice, we used the Bulirsch-Stoer method p2] to integrate the 
coupled, nonlinear, differential equations from j/i and Bode's five point method 
[p2| to evaluate integrals. Relative errors were controlled to within 10^^^. 
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Figure 1: The extremum solution curves at different values of F = 27//? (from 
the top to the bottom F = —0.9, —0.5, 0, 2, 5.) Note that each curve has an upper 
bound (3^ (F) , and no extrema can be found with f3 > P^ (F). Also note that 
over a certain range of /?, multiple solutions are present. However from those 
multiple solutions, only those are stable which have the highest value of ^0 (See 
Figure 2). The asymptotic cases are ^0 ^ — 00 (/? ^ /^o ^^'^ Pv ('") ^ 1/^^) 
and \l/o — >■ 00 which is the high temperature limit. 

6 Numerical results 

The numerical solutions of (p4) are presented in Figure 1 for the case of the 
spherically symmetric point mass system. The numerical results for the other 
models are qualitatively very similar, so we will use the point mass system to 
demonstrate their general features. As we can see, the solutions are bounded 
and, as a comparison, the results for the familiar isothermal sphere model (F = 
0)1 21 |8| are presented as well. Examination of Figure |l| clearly demonstrates 



the existence of an upper bound, j3^ (F) , which means that below a critical 
temperature there's no extremum solution at fixed F. In fact, this can be proved 
rigorously from the differential equations. ||2^ We can also see from the graph 
that J3 = -^ < and that, when /3 ^ /3q = hm /3 (^0) j the number of 

solutions goes to infinity. 
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Figure 2: The extremum solutions for F = —0.74. For f3 — 2.7, three extremum 
solutions are present, but only the one with the largest F is stable. For compar- 
ison, also see Figure 8. 

In Figure 0, we show plots of three solutions for a given temperature (we 
don't know so far which of them, if any, are stable), and in Figure we plot 
the volume density profiles of these solutions which have been normalized to the 
central density. As we can see, solutions represented by smaller ^o a-re more 
and more concentrated at the center. As a comparison, we also give the density 
profiles of the well known isothermal sphere (F = ). There's some difference 
in the shape of the density profiles in the case where F ^ but, in general, the 
volume density profile is singular as ^o ^ ^oo, and has the asymptotic solution 
Pv = T^ • Note that the linear density is p = 1 in this asymptotic case. 
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Figure 3: The relative volume density profiles of the three extremum solutions 
in the case of the isothermal sphere (T = 0, (3 — 2). Only the least concentrated 
solution, which has the highest value of '^q, is stable. The lower the value of ^o, 
the more condensed and unstable the solutions become. The unstable solutions 
are characterized by a pronounced core-halo structure. 



In order to see the difference between density profiles of the isothermal sphere 
( r = 0) and the others, we plot the high temperature solutions (Figure ^. 
These are the ^o ~^ co asymptotic solutions where /3 -^ 0. If F > the relative 
volume density profiles curve down, but when F < the profiles curve up, 
indicating that if we have an P reservoir, at higher radius the density profile 
should change from the homogeneous density profile. 
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Figure 4: In order to show the effect that a nonzero value of 7 has on the 
density, we plot the relative volume mass density profiles of the locally stable 
high-temperature asymptotic solutions (^0 — > 00) at different values of F (from 
the top to the bottom, F = —0.5,0,2,5). While the isothermal sphere has a 
constant density in this limit, nevertheless the models with F 7^ have different 
behavior due to the non-zero 7. 



In order to show how the L2 constraint affects the shape of the density 
profiles, in Figure. Q we plot the relative volume density profiles for a high 
temperature. As we can see in the figure, when 7^0 the density profiles arc no 
longer uniform and, depending on the sign of 7, the density is cither increasing 
(F < 0) or decreasing (F > 0) while in the standard case (F — 0) in the limit of 
high temperature, the density is uniform. We can understand this behaviour 
if we recognize that in the limit /3 — > gravity can be neglected. Consider 
|lO| when 7 7^ 0. Since the kinetic energy contribution is still Maxwellian, the 
probalility of finding a particle with large / is smaller when F > (7 > 0) , which 
means fewer particles will occupy larger radii. From the physical point of view, 
for a relatively large L2, more particles should concentrate at larger radii in 
order to maintain the large value, while for relatively small L2, fewer particles 
should settle at large radii in order to balance the centrifugal forces. For the 
case where F < the situation is the opposite, and the density profile should 
increase with increasing radius. Of course, while we cannot use this argument 
for finite temperatures, the origin of the difference in the density profiles when 
F 7^ is due to this effect . At finite temperature the tendency persists but the 
behavior is not guaranteed. 
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Figure 5: The f] vs E plot of the extremum solutions when L2 — 0.25 {E — 
L2 ensemble). As seen, there are no solutions below a critical value of the 
energy Ec{L2), and multiple solutions are present for a certain range of E. 
Based on the stability investigations, only the upper envelope of the complete 
spiraling curve F ^ represents the stable solutions. Other pieces are unstable 
(more eigenvalues become negative) as we go to the center of the spiral. For 
comparison, see Figure 7. 

In Figures g and || the stability properties are presented according to the 
previously discussed Poincare's linear series of equilibria with fixed L2 and E 
in the E ~ L2 ensemble. As we can see, the gravothermal catastrophe holds for 
the E ~ L2 ensemble as well. In Figure ^ below a critical energy there's no 
extremum solution (here the radius is fixed, not the energy as iny). Also, at 
the same point, (3^ (Ec) = 00 jumps to /3f, {Ec) — —00. ( We follow the spiral in 
the counter-clockwise direction.) 



18 



Figure 6: The 7 vs L2 plot of the extremum solutions when E — — 0.3.(-E — 
L2 ensemble). Above a critical value of £2, there are no extrema. In this 
case, the Gravothermal catastrophe is also present according to the stability 
investigations, and only the first piece of the spiraling curve which starts at 
L2 — and ends at the critical value of L2 defined above represents the stable 
solutions. 



Thus only the branch labelled a is stable, and the other branches are increas- 
ingly unstable because additional eigenvalues become negative. In the opposite 
case, when we fix E (Figure ||), we can see that above a critical value of L2 
there's no extremum solution. The corresponding functions, S {E) and S {L2) , 
are plotted in Figure |^. We see that the entropy is monotonic up to a critical 
energy in the stable region of the extremum solutions. 
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Figure 7: The entropy vs. energy curve for the case of L2 = 0.25 {E — L2 
ensemble). (Also see Figure 5.) The stable branch of the solutions in Figure 5 
can be identified as the highest entropy curve segment in the S vs E plot, which 
exactly ends at Ec- Note that, this segment is also a monotonic function of E, 
which means no phase transition is present in the system. 

In the T — 7 ensemble, the results are presented in Figures ^ and g. First 
of all, in the case of fixed 7, we have to take a close look at the solutions in 
the (— -E) — (3 graph. Only the first branch of the solutions are locally stable 

up to a critical value of (3 , say (3^ because i^Ej (/3^) = c», and there's no 

extremum for (3 > (3^. As a comparison, we selected a value of F in Figure m 
corresponding to Figure @; as we wind along the curve to a particular /3g value, 
the extremum solutions are represented by a more concentrated set of density 
profiles. Another result of the stability investigations is that in Figure m only 
those solutions are locally stable which are to the right of the largest maximum 
of the (3 (^0) curves. The corresponding -F can be seen in Figure lOl In the 
stable region, it's a monotonic function, and therefore there's no sign of a phase 
transition. The same holds in the E — L2 ensemble as can be seen by inspecting 
the entropy curves in Figure [7[ The results for fixed /? are shown in Figure y. 
The stable region of extremum solutions becomes unstable at 7^, and at 7 > 7^ 
there are no extrema. The free energy behaves similarly to the case of fixed 7 : 
it's monotonic, and there's no phase transition. 
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7 Conclusion 

The main purpose of this work is to study the equihbrium thermodynamics 
of sphericaUy symmetric self-gravitating systems in the mean field limit. We 
investigated both the spherically symmetric point mass system and shell sys- 
tems of differing dimension confined in a sphere. These systems are related to 
each other both in the form of their Hamiltonian and their equilibrium states. 
Furthermore, the three dimensional shell system has the interesting and poten- 
tially useful property that, with regard to the equilibrium density profile, it is 
equivalent to the spherically symmetric point mass system. Our description is 
more general than the standard treatment of the isothermal sphere pl|| since, 
in addition to the energy, we take into account L2, the sum of the squares of 
the individual angular momentum, which is conserved in the mean-field limit 
for spherically symmetric systems. In this new type of microcanonical descrip- 
tion {E — L2 ensemble) we evaluated the " equilibrium" one-particle probability 
density function for each type of system by finding the extrema of the entropy. 
The resulting pdf 's turned out to be similar to Eddington's anisotropic density 
function. |16| Therefore the density profiles obtained here also differ from the 
isothermal sphere which, in our formulation, is the special case 7 = 0. Near 
the system center, the density profiles are similar to those of the isothermal 
sphere. However, as the outer boundary is approached, depending on the value 
of r, deviations can become large, increasing or decreasing depending on the 
sign of r. The physics behind this behavior is simple: if the system is spun up 
corresponding to negative F, the outer density increases. If, on the other hand, 
the radial kinetic energy dominates the rotational energy, F > and the outer 
density decreases. 

In addition to the microcanonical ensemble, the Lagrange multiplier 7, which 
arose from the constraint on L27 yielded a new type of canonical ensemble (T— 7) 
which corresponds to the system being embedded in a heat bath at temperature 
T and an P reservoir at 7. But this analogy is only correct if the system is, at the 
least, in a local equilibrium state, and this issue demonstrates the importance of 
checking the stability properties of the extremum solutions. The method known 
as Poincare's linear series of equilibria proved adequate to analyze the stability 
of the extremum solutions in both ensembles with a little extra effort. Using 
it, we showed that only certain types of solutions are locally stable, and others 
are saddle points. In other words, the gravothermal catastrophe is also present 
both in the E — L2 and T — 7 ensembles, which means that there is no phase 
transition in our spherical mean- field models, although one was expected to be 
present because of the L2 constraint. 

From their description of the gravothermal catastrophe it is easy to imagine 
that Lyndcn-Bell and Wood had in mind a dynamical process in which mass 
was transferred from the halo to a concentrated central core. However, their ap- 
proach was confined to a comparison of stationary states. Hints of collapse have 
been seen in some N particle simulations. p4| In future work we plan to study 
the complete dynamics of the collapse, both analytically and using dynamical 
simulation of N particle and N shell systems. The shell systems should prove 
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Figure 8: The {—E) vs (3 plot of the extremum solutions in the T— 7 ensemble at 
7 = — 1 . There are no extrema above the critical (3^ (7) , and only the first branch 
of the extrema are stable. That piece of the spiraling curve which represents 
the stable solutions starts at /? = 2 and ends at (3^ (Gravothermal catastrophe 
in the T — 7 ensemble) From Figure 2, we can clearly see that only the least 
condensed density profiles represent the stable solutions. For comparison, see 
Figure 10. 
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Figure 9: The (— ^2) vs 7 plot of the extremum solutions at /3 = 2 (T — 7 
ensemble). There are no extrema above the critical value, 7^, (/3) , and only the 
first branch of the extrema are stable. The curve segment which represents 
the stable solutions starts at 7 = — 1 and ends at 7^. In other words, the 
Gravothermal catastrophe is present. 
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Figure 10: (—F) vs f3 plot at fixed 7 = — 1, where F is the thermodynamical 
potential in the T — 7 ensemble. From Figure 8, we can clearly identify which 
piece of the curve represents the stable solutions, which is the first piece, starting 
at /3 = 2 and ending at /3^. There is no sign of a phase transition because in 
the stable region —F {(3) is monotonic. Note that the stable solutions have the 
minimal F. 
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especially useful since they avoid the complication arising from the formation 
of tight binaries B. 
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